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ABSTRACT 


The  use  of  bivariate  gamma  distributions  in  simulation 
modeling  is  considered.  A family  of  algorithms  is  given,  any 
member  of  which  can  be  used  to  generate  random  bivariate 
vectors  having  any  gamma  marginal  distributions  and  any 
theoretically  possible  correlation,  both  positive  and  negative. 
Computational  considerations  are  discussed,  including  the  selection 
of  parameters  to  provide  regression  curves  consistent  with  the 
modeler's  needs.  A modification  which  is  easier  to  implement 
and  provides  most,  but  not  all,  correlations  is  also  given. 

Finally  the  use  of  these  algorithms  to  generate  first  order 
autoregressive  time  series  is  discussed. 

<\ 
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Many  random  phenomena  may  be  modeled  as  gamma  random  variables 
and  many  random  phenomena  are  correlated.  For  example,  time  between 
"events"  in  the  analysis  of  nulcear  power  plant  safety,  population 
characteristics  such  as  income  and  age,  and  species'  lifetimes  are  all 
nonnegative  random  variables  which  may  be  modeled  by  the  gamma  distri- 
bution in  some  cases.  In  addition,  these  variables  are  often  subject 
to  the  same  exogenous  variables,  causing  them  to  be  either  positively 
| or  negatively  correlated.  For  example,  two  power  plants  may  be  affected 

by  the  same  earthquake  causing  positive  correlation,  age  and  income  may 
be  correlated,  and  a higher  than  average  rainfall  may  be  helpful  to  one 
species  and  harmful  to  another  species,  causing  negative  correlation  in 
their  lifetime. 

While  it  would  then  appear  that  bivariate  or  multivariate  gamma  models 
should  appear  frequently  in  the  literature,  this  does  not  seem  to  be  the 
case.  Instead,  the  multivariate  normal  distribution  is  commonly  used 
when  correlated  random  variables  are  appropriate,  since  multivariate 
normal  properties  are  well  known  and  random  vectors  can  be  generated 
relatively  easily  (Scheuer  and  Stoller  [34],  Deik  [9]  and  Schmeiser  and 
Ali  [35]). 

In  the  unidimensional  case,  the  gamma  distribution  is  a generaliza- 
tion of  the  normal  in  that  the  normal  distribution  is  a limiting  distri- 
bution of  the  gamma  as  the  shape  parameter  goes  to  infinity.  Thus  in  the 
multidimensional  case,  a multivariate  gamma  model  is  more  general  than 
the  multivariate  normal  model,  thereby  allowing  more  general  modeling 
ability. 

In  this  paper  a bivariate  gamma  distribution  is  developed  in  the 
form  of  a family  of  algorithms  capable  of  generating  random  vectors 
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possessing  gamma  marginal  distributions  and  specified  correlation 
coefficient  p.  In  particular  the  marginal  distributions  have  den- 
sity functions 

V1 

f±(x)  - ((x/e±)  1 exp(-x/3±)  / (g±  r(a1)) 

for  x > 0,  a±  > 0,  3±  > 0,  i - 1,  2 

The  correlation  p does  not  depend  upon  the  scaling  parameters  8^ 
and  $2»  so  we  assume  8^  = $2  = ^ in  t*ie  rest  this  paper.  The 
desired  scaling  can  be  achieved  by  multiplying  the  generated  variates 
by  their  respective  scaling  parameters,  as  is  done  in  the  algorithm 
in  the  Appendix. 

The  results  in  this  paper  differ  from  previous  work,  discussed 
in  Section  1,  in  that  the  algorithms  developed  here  do  not  require 
integer  shape  parameters  a^,  nor  equal  shape  parameters,  and  they 
provide  variates  having  any  theoretically  possible  correlation  p. 

As  seen  in  later  Section  1,  most  existing  methods  allow  only  positive 
correlations,  and  then  not  all  possible  positive  correlations. 

Thus  the  range  of  possible  correlations  for  given  shape  parameters 
and  is  of  interest.  For  the  bivariate  normal  p £ [-1,  1] . How- 
ever for  gamma  marginal  distributions,  not  all  correlations  are  consistent 
with  particular  shape  parameter  values.  Figure  1 shows  the  obtainable 
correlations  as  a function  of  a2»  given  * 1 and  5.  Correlations 
inconsistent  with  the  specified  shape  parameters  are  shaded.  Note  that 
only  when  * a 2 is  it  possible  for  p ■ 1.  Likewise  p = -1  is  not 
possible  except  in  the  limit  as  or  a2  tend  to  infinity.  The  dashed 
lines  are  discussed  in  Section  1. 
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Figure  1 About  Here 


The  maximum  and  minimum  possible  correlations,  given  in  Moran 


F'^l-F  (X  )),  re- 
a2  ax  1 


[31],  occur  when  X ■ F 1 (F  (X  ))  and  X. 

a ^*1  * * 

spectively,  where  F^Cx)  and  F^Cu)  are  the  cumulative  distribution 
function  (cdf)  and  inverse  cdf,  respectively,  of  the  gamma  distri- 
bution with  shape  parameter  a.  These  relationships  are  symmetric 
and  can  alternatively  be  expressed  by  interchanging  all  subscripts. 
In  Figure  2,  these  relationships  are  illustrated  for  a specific 
value  = x^,  with  = x2  being  the  corresponding  value  when 
correlation  is  maximized  and  X2  = x2  being  the  corresponding  value 
when  the  correlation  is  minimized. 


Figure  2 About  Here 


1.  EXISTING  METHODS  FOR  BIVARIATE  GAMMA  GENERATION 


There  are  several  existing  methods  for  generating  random  bi- 
variate gamma  vectors,  although  each  has  a limitation  which  makes 
it  less  general  than  the  algorithms  developed  in  Section  2. 

Trivariate  reduction,  developed  by  Arnold  12], -can  be  used 

to  generate  variates  having  positive  correlation  with  p < min  (a^, 

1/2 

0t2)/ (cijC^)  . Letting  "gamma  (a,  g)"  denote  the  gamma  distribution 
with  parameters  a and  g,  the  trivariate  reduction  algorithm  proceeds 
as  follows  for  given  parameter  values  a^,  o^,  end  p. 

ALGORITHM  GTVR 

1/2 

1.  Generate  Y1  ^ gamma(a^  - p(a^a2>  , 1) 

1/2 

2.  Generate  Y2  'v  gamma  (a.,  - pto^o^)  , 1) 

1/2 

3.  Generate  Y^  ^ gamma (pCa^o^)  , 1) 

4.  X2  «-  Yx  + Y3 

5 . X2  «-  Y2  + Y3 

As  noted  earlier,  X^  and  X2  can  be  multiplied  by  g^  and  g2,  respectively, 
to  obtain  any  desired  scaling.  In  algorithm  GTVR,  Y3  is  a common  compo- 
nent of  both  X^  and  X2>  thus  inducing  positive  correlation.  The  dashed 
lines  in  Figure  1 indicate  the  boundaries  of  correlation  obtainable  with 
trivariate  reduction.  When  a positive,  but  not  extreme,  correlation  is 
needed,  trivariate  reduction  is  an  excellent  algorithm.  Fishman  [11] 
suggests  this  method  for  generating  correlated  gamma  variates. 

Ronning  [33]  generalized  trivariate  reduction,  although  no  specific 
algorithms  are  given.  Rather  a general  framework  is  discussed,  from  which 
it  appears  possible  to  create  algorithms  which  would  be  able  to  obtain  any 
positive  correlation  p.  In  addition,  he  considered  the  n dimensional  case. 


4 


*1 


Composition,  or  probabability  mixing,  provides  a very  general 
algorithm  for  generating  random  bivariate  vectors,  in  that  all  possi- 
ble parameter  values  may  be  obtained,  including  negative  correlations. 
The  algorithm  proceeds  as  follows: 


ALGORITHM  GCOMP 

1.  Generate  U 'v*  11(0,1) 

2.  If  U > p/C,  go  to  step  5. 


3.  Generate  'v*  gamma(a^,  1) 

4.  Generate  ^ gamma (a2»  1) 

5.  Generate  V ^ U(0,1)  (or  V 

6.  X,  *-  F-1(V) 

1 

7.  If  p < 0,  V «-  1 - V 

8.  X «-  F‘1(V) 

2 a2 

9.  Deliver  (X^  X2> 


and  go  to  step  9. 
■ (U*C  -p)/(C  -p)) 


Here  C is  the  maximum  or  minimum  possible  correlation  given  and  a^. 

C is  calculated  using  equation  (3)  in  Section  2 with  H^  M a^,  Hj  * ot 
and  y ■ 0.  The  implementation  of  algorithm  GCOMP  could  proceed  in  many 
ways,  with  the  given  logic  being  an  example.  The  underlying  idea  is  to 
generate  an  independent  vector  (X^,  X2>  with  probability  p/C  and  to 
generate  a vector  (X.,  X2)  having  correlation  C otherwise. 

The  most  severe  problem  with  the  composition  approach  is  the  sta- 
tistical properties  of  the  generated  variates.  Figure  3 is  a scatter 
plot  of  1000  random  generated  vectors  with  ot^  * 10,  a2  ■ 10,  and  p =“.5. 
The  marginal  distributions  and  the  correlation  are  as  desired.  However, 
the  points  fall  in  a pattern  that  is  inappropriate  in  most  applications. 
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The  pattern  is  that  of  many  independent  variates  superimposed  over 

many  other  variates  lying  on  the  curve  X„  = F ^ (l-F  (X.)). 

^ (*2  i 

Examination  of  Figure  3 thus  shows  that  in  addition  to  obtaining  the 
correct  marginal  distributions  and  correlation,  some  thought  must  be 
given  to  the  pattern  of  vectors  obtained. 


Figure  3 About  Here 


A third  approach,  advocated  by  Kottas  and  Lau  [25],  attacks  the 
pattern  directly,  by  modeling  the  regression  line  E(X2[x^)  and  the 
variance  V(X2|x^)  explicitly.  Their  work,  done  in  a much  more  general 
context,  does  not  provide  for  specified  marginal  distributions  or  spe- 
cified correlations.  Thus,  to  a large  extent,  the  method  chosen  depends 
upon  the  data  available  to  the  modeler.  When  the  structure  of  the 
dependency  between  X^  and  X2  is  understood,  using  this  understanding 
is  an  excellent  approach.  In  related  work,  Hull  [17]  approximates  the 
desired  correlation  while  considering  the  regression  line. 

Moran  [31],  Gunst  and  Webster  [16],  Johnson  and  Kotz  [22]  and  Kibble  [23] 
discuss  and  reference  several  other  multivariate  gamma  distributions, 
but  all  have  restrictions  such  as  = ot^  , 2a  ^ integer,  or  p > 0. 

While  none  were  developed  for  use  in  simulation,  they  could  be  useful 
in  some  cases. 

There  are  several  papers  which  provide  methods  which  could  be  used 
to  approximate  a bivariate  gamma  distribution  with  specified  correlation 
p > 0.  Lee  [26]  discusses  a multivariate  Wei’jull  distribution.  Takahasi 
[39]  and  Durling,  Owen,  and  Drane  [10]  discuss  multivariate  Burr 
distributions  with  p > 0.  Mihram  and  Stacy  [36]  discuss  a warning-time/ 
failure-time  bivariate  distribution  with  beta  and  generalized  four-parameter 
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gamma  marginal  distributions.  Johnson  [21]  considers  translations  of 
bivariate  normal  vectors.  Johnson  and  Ramberg  [19],  developed  a bi- 
variate uniform  distribution. 

For  p < 0,  few  methods  have  been  developed.  Gargano  and  Tenenbein 
[13]  and  Johnson  and  Tenenbein  [20]  have  developed  bivariate  uniform 
distributions  which  can  be  used  tc  generate  variates  having  negative 
correlations.  Using  rank  correlation,  bivariate  gamma  distributions 
can  be  generated  by  transforming  the  bivariate  uniform  variates  to  gamma 
variates  via  an  inverse  cdf  transformation. 

In  the  related  area  of  autocorrelated  sequences,  two  papers  suggest 
methods  for  p < 0:  Polge,  Holliday,  and  Bhagavan  [32]  suggest  a "corre- 
lation transfer"  method,  which  provide  ■:<  approximate  correlation  and 
marginal  distributions;  and  Gaver,  Lavenberg,  and  Price  [14]  propose 
generating  correlated  variates  with  negative  correlation  via  a Markov 
chain  model. 


2.  A FAMILY  OF  ALGORITHMS 


Developed  in  this  section  is  a family  of  algorithms,  any  member 
of  which  can  produce  bivariate  gamma  vectors  having  any  parameters 
c^,  a2  and  p. 

Proposition  1.  If  Z,  W^,  and  are  independent  gamma  random  variables 
with  shape  parameters  y,  6^,  and  62>  respectively;  and  U is  an  independ- 
ent U(0,1)  random  variable;  and  either  VeUorV*l-0;  then 

X1  ■ FhJ<0)  + Z + H1 

and  (1) 

X2  . F’V)  + Z + W2 

are  each  gamma  random  variables,  with  shape  parameters  + y + 6^ 

and  a2  * H2  + y + 62>  respectively. 

Proof:  The  result  follows  immediately  from  the  reproducibility  property 
of  the  gamma  distribution  and  from  noting  that  F ^(U)  and  F 1(1-U)  are 
each  random  variables  having  cdf  F. 

Proposition  2.  Corr^,  X2>  «[E(F"1(U)F‘1(V))-H1H2  + yJ/tejC^)172 
Proof:  The  proof  is  straightforward  algebra. 

Propositions  1 and  2 say  that  equations  (1)  can  be  used  to  generate 
(X^,  X2)  having  the  desired  marginal  distributions  and  correlation. 

Given  and  a2>  any  theoretically  possible  correlation  p may 
be  obtained.  Examination  of  some  important  limiting  cases  is  helpful 
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to  understand  the  robustness  of  the  algorithms: 


1.  Let  y = 6^  = 6^  = 0.  Then  the  Corr(X^,  X^)  is  maximized  if 
V * U and  minimized  if  V = 1 - U,  as  discussed  earlier. 

2.  Let  1 0 or  = 0,  and  y «=  0.  Then  X^  and  X^  are 
independent  random  variables. 

Letting  the  parameters  range  between  these  two  extreme  cases  yields 
all  possible  desired  correlations. 

The  remaining  problem  is  to  select  values  of  the  five  parameters 
Hl’  H2’  Y*  *^1  an<^  *^2  t0  °^tain  desired  marginal  distributions  and 
correlation.  The  following  conditions  must  be  satisfied: 


H1  + Y + 61  = “l 
H2  + y + 62  = a2 

[E(Fy1(U)F~^(V))-(H1H2)  + y]  / (a^)172  ■=  p 
Hlf  H2,  y,  61>  62  > 0. 


(2) 


Since  we  are  using  five  variables  to  satisfy  three  equality  conditions, 
finding  a set  of  parameter  values  corresponds  to  finding  a feasible 
solution,  rather  than  an  optimal  solution,  to  a nonlinear  programming 
problem.  Here  6^  and  $2  are  slack  variables. 

An  efficient  solution  procedure  for  determining  parameter  values 
is  important,  since  substantial  computation  is  required  to  determine 
whether  or  not  conditions  (2)  are  satisfied  for  given  parameter  values. 
Most  of  the  computation  is  involved  in  calculating 

Corr(X1,  X2)  = [E(F~V)  f'^V))-^^)  + y]  /(a^)172  (3) 

since  the  expected  value  must  be  calculated  numerical.' y using  any  one 
of  the  following  three  integrals: 
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(4) 


F "V)  F ~\u)  du 

0 H2  h1 


F *(F  fx))  x Hlexp(-x)  dx  / T(H.)  (5) 

0 H2  Hx  1 


F u^(FH(-ln  y))  (-In  y)  Hldy  / r(H.)  (6) 

0 ’ 1 1 


if  p > 0.  If  p < 0,  then  replace  F F(u)  with  F F(l-u)  in  equation  (4) > 

H1  H1 

replace  F (x)  with  1-F  (x)  in  equation  (5) , and  replace  Fu  (-In  y)  with 

H1  H1  H1 
1-F  (-In  y)  in  equation  (6).  Davis  and  Rabinowitz  [8]  and  other 
H1 

numerical  integration  textbooks  discuss  various  methods  for  evaluting 
these  integrals.  We  seemed  to  have  best  results,  in  term."  of  a sub- 


jective trade  off  between  speed  and  accuracy,  using  a 24  point  Gaussian 


method  on  integral  (4).  Integrals  (5)  and  (6)  have  the  advantage  of 


requiring  only  one  inverse  cdf  evaluation,  which  is  advantageous  since 
F '*'(u)  is  iteratively  calculated  from  F(x)  in  almost  all  published 


methods.  (See,  e.g.,  Rest  and  Roberts  [5]). 


One  way  to  select  the  parameter  values  is  to  choose  the  feasible 
values  satisfying  conditions  (2)  and  to  optimize  by  making  the  curves 
of  regression  E(X^|x2)  and  E(X2|x^)  behave  as  desired.  Proposition  3 
gives  these  curves  as  a function  of  the  parameters. 


Proposition  3. 

fl 


e(x2|x1)  = 


If  X^  and  X2  are  defined  as  in  equations  (1)  and  p > 0, 

-i  V1  “r1 

FH  (FH  ^X1  u (l“u)  76(1^, o^)  du  + Yx1/a1  + 62 


and  E(X2|x2)  is  obtained  symmetrically  by  interchanging  all  subscripts 
"1"  and  "2".  If  p < 0,  replace  Fu  (x.u)  by  1-F„  (x,u). 
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Proof:  The  results  follow  directly  given  that  the  distribution  of 
F^(U)  given  x is  beta  with  parameters  (H^,  a^)  over  the  range  [0,  x] 
and  the  distribution  of  Z given  x is  beta  with  parameters  (y,  a^)  over 
the  same  range. 

An  important  special  case  arises  when  and  p > 0,  since  then 

e(x1|x2)  = x2(h1  + y)/«2  + 51 

and 

E(X2|x1)  = x1(H2  + y)/a1  + 62. 

Note  that  these  lines  are  not  identical,  nor  are  they  even  parallel, 

1/2 

although  they  may  be  made  parallel  by  setting  y = (a^a2)  - H^. 

The  ability  to  obtain  nonlinear  lines  of  regression  is  an  important 
advantage  of  this  general  family  of  algorithms  compared  to  the  more 
limited  trivariate  reduction  algorithm,  in  which  always  = H2  = 0. 

The  generation  of  the  bivariate  vector  (X^,  X2>  can  be  performed 
directly  from  Equations  (1),  with  Z,  W^,  and  V?2  being  generated  by  any 
gamma  generator.  Theoretically  exact  gamma  generators  are  given  in 
Ahrens  and  Dieter  [1],  Atkinson  [3],  Atkinson  and  Pearce  [A],  Cheng  [6], 
Fishman  [11,  12],  Johnk  [18],  Greenwood  [15],  Kinderman  and  Monahan  [24], 

Marsaglia  [27],  McGrath  and  Irving  [28],  Schmeiser  and  Lai  [36],  Tadi- 
kamalla  [37,  38],  C.  S.  Wallace  [40],  N.  D.  Wallace  [41],  and  Whittaker  [42]. 
Speed  is  not  a factor  here  since  the  evaluation  of  F ^(U)  requires 
orders  of  magnitude  more  time  than  the  fastest  of  these  algorithms. 

Care  need  be  taken  only  to  not  use  an  algorithm  with  time  requirements 
which  are  unbounded  for  some  parameter  values. 

Considerable  execution  time  can  be  saved  by  implementing  the  algorithm 
less  directly.  Rather  than  performing  two  slow  evaluations  of  F 1 , the 
following  requires  only  one  F ^ evaluation: 
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Generate  ^ gamma (H^) 

0 * FHl«l> 

If  P < 0,  u v 1 - u. 

Generate  Z ■v  gamma  (y,  1) 
Generate  'v-  gamma (6^,  1) 

Generate  ^ gam.ua  (6  ^ , 1) 

«-  xx  + Z + w 

X2  "■  FH1(U)  + 2 + w2 


3.  SELECTION  OF  A PARTICULAR  ALGORITHM 


Several  criteria  can  be  identified  for  selecting  a particular 
algorithm  from  the  family  developed  in  Section  2: 

1.  Set-up  time  required. 

2.  Marginal  execution  time  for  each  vector  (X^,  X^)  generated. 

3.  Program  bulk  (related  to  lines  of  code  and  support  routines 
needed) . 

and 

4.  Statistical  properties  of  the  generated  variates. 

This  fourth  criterion  is  probably  the  most  important  since  the  first 
three  criteria  are  very  similar  for  all  algorithms  in  the  family.  Thus 
we  chose  one  algorithm  for  implementation  based  on  criterion  4,  as  mea- 
sured subjectively  by  viewing  scatter  plots  for  many  values  of  a^, 
and  p.  The  algorithm  given  in  detail  in  the  appendix  seemed  to  behave 
as  "expected".  While  none  of  the  algorithms  were  "bad",  some  tended  to 
produce  a rather  sharp  boundary  on  the  generated  values.  The  algorithm 
selected  tended  to  produce  values  which  tapered  off  gradually,  which  to 
the  authors  appeared  to  be  subjectively  better.  Figure  4 is  a scatter 
plot  produced  by  the  selected  algorithm. 

Figure  4 About  Here 

The  algorithm  implemented  in  the  Appendix  corresponds  to  y = 62  ~ 0 
and  to  assuming  that  a2’  The  second  assumption  is  without  loss  of 

generality  and  is  implemented  in  the  algorithm,  making  it  transparent 
to  the  user. 
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0 1 H 1 °1 

[E(F"1(U)F~1(V))  - H^2J/(a1a2)1/2  = p 

Determining  the  value  of  H satisfying  conditions  (8)  is  accomplished 
by  bounding  H and  then  using  the  modified  regula  falsi  root  finding  algo- 
rithm. (See,  for  example,  Conte  and  deBoor  [7].)  The  bounds  used  were 
good  enough  that  never  more  than  four  iterations  are  needed  to  obtain  p 
within  .0001.  On  the  average  two  or  three  iterations  are  required.  The 
number  of  integral  evaluations  is  two  more  than  the  number  of  iterations. 

The  bounds  are  given  in  the  following  propositions.  cmin  and  C ^ 
are  the  minimum  and  maximum  correlation  theoretically  possible,  calcula- 
ted from  equation  (3)  with  ■ a^,  H2  ■ a2>  and  y » 0. 

Proposition  4.  If  p « C . , then  H * a, . 

mm  1 

2 

Proposition  5.  If  C , < p < 0,  then  a, (p/C  . ) < H < a, . 

min  l min  1 

Proposition  6.  If  p * 0,  then  H K 0. 

1/2  2 

Proposition  7.  If  0 < p < (a^/o2>  , then  a^p  < H < a^. 

1/2 

Proposition  8.  If  p * (o^/a^)  , then  H - a2« 
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[ 

1/2  3 

Proposition  9.  If  (a./a  ) < p < C , then  a.  < H < a. (p/C  ) . 

2 1 max  2 1 max 

Proposition  10.  If  p ■ C , then  H • a, . 

max  1 

The  proofs  follow  directly  from  Figure  5. 
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4.  A RELATED  METHOD 

In  this  section  a method  is  given  which  will  generate  negatively 
correlated  gamma  variates  with  closed  form  determination  of  the  appro- 
priate parameter  values.  It  will  not,  however,  provide  all  possible 
negative  correlations. 

Similar  to  equations  (1),  consider 


Xx  - Y1  + Z + Wx 


Y2  + Z + W2 


(9) 


where  Z,  W^,  and  have  independent  gamma  distributions  with  shape 
parameters  y,  6^  and  62;  and  and  are  Erlang  k distributions 
obtained  as 


k 

Y - -In  (tt  U ) 
i-1  1 

and  . 

k 

Y2  - -In  (tt  (1-U  )) 
i*l 


where  the  U^'s  are  independent  U(0,1)  random  variables.  The  correlation 
of  Y^  and  is  _-64493  for  all  values  of  k.  Thus  Corr(X^,  Xj)  ■ 
(-.64493  k + y)  / (a^)172. 

Given  a^,  a^,  and  the  correlation  of  X^  and  X2»  P,  the  appropriate 
parameters  can  be  determined  as  follows: 


k INT  (p(a1a2)1/2  / (-.64493)) 

Y «■  -.64493  k + p(aLa2)1/2 

61  * “l  ” k " Y 

62  ^ a2  " k " Y 

If  k > 1;  and  Y.  ^ and  62  are  nonnegative;  the  desired  corre^.ion  is 
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achieved  via  equations  (9) . Otherwise  the  methods  of  the  previous  two 
sections  can  be  used. 


5.  EXTENSION  TO  TIME  SERIES  MODELS 

Consider  the  generation  of  a random  time  series  X^, 
where  each  X^  has  a gamma (a,  B)  distribution  and  Corr(X^  * Xi+1> 
An  autoregressive  time  series  having  these  properties  can  be 
generated  using  the  algorithm  of  Section  3 as  follows: 


ALGORITHM  GTS 


1 . Generate 


^ gamma(a,  1)  • X^  ■*-  B 


2 . For  i - 2 , 3 , 

Generate  ^ gamma (a-H,  1) 

v * VW 

If  p < 0,  V 1 - V. 

YI  * r H1(V)  + W1 


*£*  6 Y£ 


where  H is  calculated  to  satisfy  equation  (8)  with  - «2  - a* 

The  higher  order  correlations  CorrCX^,  X^+^)  will  tend  to  zero 
as  k becomes  large,  but  the  calculation  of  the  specific  values  is 
not  closed  form.  Extension  of  the  algorithms  of  Section  2 follows 
in  a similar  manner. 
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APPENDIX 


■ 

This  appendix  gives  logic  to  determine  parameter  values  and  to 
generate  random  bivariate  vectors  (X^,  X2)  using  the  algorithm  developed 
in  Section  3.  Parameters  which  are  provided  are  cx^,  a^,  B2»  P and 

EPS.  EPS  is  the  error  allowed  in  p when  determining  H.  Corr(H,  a^, 
a^t  p)  is  the  correlation  of  X^  and  X^  given  the  specified  parameter 

) 

values  and  is  obtained  by  equation  (3)  which  involves  numerical  inte- 
gration of  (A),  (5),  or  (6)  in  Section  2 with  ■ H,  ■ a^,  and 

Y * 0.  Care  should  be  taken  to  use  an  integration  routine  which  will 

1/2 

evaluate  (A),  (5),  or  (6)  with  error  less  than  EPStejC^) 

ALGORITHM 

Initialization  (executed  each  time  a^,  a2  or  p changes  value) 

1.  (If  p is  close  enough  to  zero,  go  generate  independent  variates.) 

If  | p | < EPS,  H 0 and  go  to  step  10. 

2.  (Ensure  that  o^.) 

IND  0.  If  < a^,  IND  ■*-  1,  T a^,  a2»  a2  *■  T . 

3.  C ■*-  Corr(a1,  o^,  a2,  p) 

If  | C | < |p|,  then  p is  not  theoretically  consistent  with  and  a2* 
Return  with  an  error  code. 

A.  (Bound  H between  A and  B with  associated  correlation  errors  F and  G 
in  preparation  for  the  modified  regula  falsi  search  in  step  6.) 

2 

If  p < 0,  A *■  a1(p/C)  , B «-  o^,  G«-p-C,  F ■*-  p - Corr(A,  c^,  a2,  p) 

1/2  2 1/2 
If  0 < p < (aj/c^)  , A •*-  a^p  , B «-  a2,  G p - 

F ■*-  P - Corr (A , a.,  a2,  p) 

1/2  2 

If  (o^/c^)  < p,  A *-  a2,  B «-  Cp/C)  , G p - Corr(B,  alf  a , p) 

F «■  p - (a2/a1)1/2 
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5. 


FO  «-  F 

6.  (Use  the  modified  regula  falsi  method  to  determine  H in  steps  6-9.) 

H «-  (G*A  - F*B)  / (G  - F) 

FN  «-  p - Corr (H,  o^,  a2,  p) 

PROD  «-  FO*FN 
FO  «-  FN 

If  F*FN  < 0,  go  to  step  8.  Otherwise  go  to  step  7. 

7.  A * H,  F ♦ FN.  If  PROD  > 0,  G ■*-  G/2.  Go  to  step  9. 

8.  B «-  H,  G FN.  If  PROD  > 0,  F F/2. 

9.  If  ABS(FN)  > EPS,  go  to  step  6.  Otherwise  H (G*A  - F*B)  / (G  - F) . 

Generation  (executed  each  time  a new  vector  is  needed) 


10.  If  H > 0 and  H < a^,  go  to  step  11. 

Otherwise  generate  gamma  (o^) , U Fq  (X^) . 


„-l. 


If  p < 0,  U * 1 - U.  Xx  «-  0.  If  H > 0,  X1  F *(U).  Go  to  step  12. 


11.  Generate  Xj*'*  gamma(H).  U FjjCXj)*  If  p < 0,  U *■  1 - U. 

X2  - F'*<U) . 

12.  If  H < a^,  generate  W gamma(a^-H)  and  X^  +■  X^  + W. 

Xx  «-  X^  X2  «-  82*  X2.  If  IND  - 1,  T Xx,  X1  *■  X2,  and  «-  T. 

Deliver  (X  , X2). 

A FORTRAN  implementation  is  available  upon  request  from  the  authors. 
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Figure.  2. 


The  relationship  between  two  random  variables 
having  minimum  and  maximum  correlation. 
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Figure  3.  A scatter  plot  of  bivariate  gamma  points  using 
the  composition  algorithm  GCOMP . 
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APPENDIX  II 


Fortran  implementation  of  the.  algorithm  developed  In  Section  3, 
Including  a driver  program  and  numerical  Integration  auhroutlne. 
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